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fj ' Abstract The altitude distribution of optical turbulence is derived from the MASS 

instrument data by solving an inverse problem. In this paper, some modifications of the 
profile restoration are described. The principal change is the introduction of the Non 
Cn , Negative Least Squares algorithm which has good regularizing properties. An averaging 

of scintillation indices was replaced with averaging of obtained solutions what leads to 
clearer physical results. It is shown that restoration with a number of turbulent layers 
as large as 14-15 can be successfully performed. 
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(Vj ' 1 Introduction 

It is well known that the efficiency of astronomical observations in optical and near IR 
fvj , range greatly depends on turbulence in the earth's atmosphere. Modern techniques for 

^ ' telescope efficiency gain require both statistically-valid long-term and near real-time 

OO I information about properties of the optical turbulence (OT) above an observatory. In 

lO ■ general, the turbulence intensity is described with help of the refraction index structure 

^P ' constant Cn- One of the instruments designed to measure this parameter is MASS 

(Multi Aperture Scintillation Sensor) [I] developed a decade ago. More than 30 such 
lO ' devices are used in different astronomical observatories and site testing campaigns [2]. 

^— ^ , The method used by the MASS 1T1|3] to determine the OT vertical profile is based on 

the simultaneous measurements of stellar scintillations in four entrance apertures {A, 
B, C, D) of different size and produces four normal and six differential [4] scintillation 
indices s — a variance of the relative fluctuations of light fluxes. The measurement 
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is performed with an integration time Te of 1 ms and an estimation of s is calculated 



from the measurements at time interval T^ (basetime) of about 1 s [3lJ5]- 

Scintillation indices s are affected by turbulence strength in the whole atmosphere 
traveled through by the light: 

f-OO 

s^= / Ci{h)Q{h)dh, (1) 

Jo 

where Cn{h) is index constant at altitude h, Q{h) is certain weighting function that 
depends on geometry of light propagation: size of receiving aperture and angular size of 
emitting object. The function Q{h) can be calculated for a given geometry and stated 
spatial spectrum of refraction index perturbation [3ll5| . 

The calculation of the vertical distribution of the OT from measured s indices is 
only possible by solving the inverse problem which is ill-conditioned for any practical set 
of MASS apertures (therefore, set of Q{h) functions). For this reason, the restoration 
algorithm of the turbulence profiles is a key part of the MASS method. The description 
of the algorithm and its verification are given in paper [5]. Instrumental and other 
factors affecting accuracy of MASS results were considered in [^Hl]. 

During development of the algorithm used in the Turbina software (hereafter — 
algorithm P), some simplifications were adopted with intuitive guesses due to lack of 
real data for careful analysis of the algorithm work. 

In recent years, intensive measurements with the MASS instrument have been 
performed by different groups at different sites. The large amount of data encouraged us 
to revise the profiles restoration technique. The revision focuses on the implementation 
of a more valid mathematical standpoint and on the verification (and correction if it 
is necessary) of the model assumptions. 

A particularly significant factor for us is the completion of a three year campaign 
at mount Maidanak [7] and a successful two year campaign at Shatdzhatmaz summit 
[H] where 2.5 m telescope of Sternberg institute should be installed. 



2 Algorithm P of Turbina software 

Formula ((T|, which comes from the theory of weak perturbation, is linear in C„(/i) 
and can be used to describe the effects of turbulence in a typical astronomical night. 
For the case of Kolmogorov's model, optical turbulence in the whole atmosphere or 
a certain altitudinal range is defined by one parameter only. Usually, Fried radius rg 
(atmospheric coherence radius) is applied, but for description of OT distribution along 
the line of sight, the turbulence intensity J is preferable 

J = [ Cl{z)dz,= Qm-\^r^^'^ or J = 1.5 ■ 10"^"* ■ r-p^^^^ if A = 500 nm (2) 

Together with J, the seeing /3 (in the conventional sense) is useful. Scaled to arcsec 
it is /3 = 2 ■ 10 ■ J ' for propagating light of wavelength A = 500 nm. 

To solve the OT restoration problem, the integral equations ([T} are written in 
discrete form for some fixed altitude grid: h = (/iq, hi ... , /i„_i): 

Ax = b, (3) 

where the vector x = (Jq, Ji . . . , Jn-l) is the solution of turbulence intensities, the 
vector b = ((sq), . . . , (sj._]^)) is a set of measured scintillation indices and A is a 



matrix fc x n of values of atmospheric weighting functions Qj{hi), calculated for given 
j-th apertures. 

Turbulence intensity Ji is represented by the integral: 

rh+Ah 

J= cl(z)dz^Ah{cl), (4) 

Jh 

where (C„) is the mean value of the structure coefficient within the i-th layer. This 
is not the only way to solve the problem, alternative methods will be discussed in 
Section [TCTl 

The standard input values that remain unchanged throughout the process are: the 
number of used indices k = 10, the number of layers n = 6, the altitude grid 0.5, 1, 2, 
4, 8, 16 km. Input data (s ) are mean values over accumulation time Ta {accumtime) 
that is about 1 minute. 

The similar characteristics of the weighting functions Qj{h) (e.g. Fig. 1 in ^) lead 
to a large condition number of the matrix A (certainly > 1000), and a poorly bound 
equation system. Without an a priori information, a proper solution of Least Squares 
Problem (LSP) can only be obtained if the input data have small relative errors which 
is encountered in the case of strong turbulence. 

The physical nature of the problem lets us apply additional restrictions: the solution 
X of the system Q must be non-negative. To solve Q with the additional restriction 
X > 0, the direct minimization of weighted sum of residual squares with the help 
of the Powell method 9 was implemented in the procedure j4irrao4J. To provide the 

1/2 

restrictions, the original variables are replaced hy v^ — J^ , i.e. the problem becomes 
nonlinear and the final solution is obtained by Uj squaring. 

The non-linearity of the problem (the function to be minimized becomes polynomial 
of 4 degree in Vi) leads to the appearance of a number of local minimums, where 
minimization process may be completed with some probability. The non-uniqueness 
of solutions, topical computational cost, and difficulties in estimating the solutions 
errors force us to find more mathematically reasonable algorithms for the OT profiles 
restoration. 



3 Non-negative least squares (NNLS) 

A solving a system of linear equations in terms of least squares with non-negative 
solutions is not new. The technique was put forward 40 years ago and named NNLS 
(Non-Negative Least Squares) [10]. The algorithm is based on the fact that the mini- 
mum of the sum of the squares R = ||Ax — b|| lies either into restricted domain or 
on the its border. To find the R minimum one cannot set a negative components of 
the non-restricted solution to zero because the main axes of the residual paraboloid do 
not coincide with the coordinates axes. It was proved that NNLS finishes in a finite 
number of iterations and usually requires only about n/2 iterations 10 . 

The NNLS technique is used in astronomical data processing, although not really 
frequently [Tnil2lll3lll4l[T5lfT6] . Goodwin [17] successfully applies NNLS for restoration 
of altitude turbulence profile from SLODAR measurements. 

Following the book [TD], we implemented NNLS and other algorithms needed to 
solve a system of linear equations related to the LSP using pure C-|— I- and boost 



^ see |http://curl. sal. msu.ru/mass/download/doc/dataproc.pdf I 




Fig. 1 Comparison of OT profiles restored by algorithm Nl (left) and algorithm P (right) 
with data taken on September 6, 2005. Turbulence intensities in 0.5, 1, 2, 4, 8 and 16-km layers 
are shown on top, free atmosphere seeing is shown at the bottom. 



(http://boost.org \ librarieqj. The LSP is solved by singular value decomposition (SVD) 



using Householder transforms and the modified QR-decomposition. 

To obtain the best unbiased solution, the source set of equations ((SJ must be 
weighted through a left multiplication by matrix W. Assuming that input errors are 
uncorrelated, the diagonal matrix W (where Wjj = (1/cr.,) ' , and (t.- is an estimator 
of variance of j— index error) is used like in the case of algorithm P. 

The replacement of the solving algorithm was the first step of the revision. This 
version is referred below as Nl. 



4 Comparison of tlie turbulence profile restoration algorithms 



In order to compare algorithms Nl and P, we used data taken with the original MASS 
device [3] at mount Maidanak between 2005 and 2007 [7]. Both algorithms were run 
with identical input parameters and data set. The data consist of more than 3 • 10 
seconds of measurements or 50 000 points. 

The principal conclusion from the comparison is that both set of solutions are 
in a good agreement, differing only in details. The results for September 6, 2005 are 
given on Fig. [T] as an example. One can see that all the characteristic features appear 
on both left and right plots. Computed with C„ profiles jHfree values are virtually 
matched. The free atmosphere seeing Pfree l5j is produced by OT above the boundary 
layer (conventionally above 1 km, here 0.5 km layer is included as well). The same 
level of agreement between the two algorithms was found on other nights regardless of 
turbulence intensity and distribution. 

The inter-comparison of the free atmosphere seeing Pfree over the whole data 
set shows that the mean-square deviation is only 0.02". The systematic difference is 
extremely insignificant, medians Pfree differ less than 0.006" between the techniques. 



LSP library is available for download http://curl.sai.msu.ru/~matwcy/lsp/ 
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Fig. 2 Normalized distributions of differences 5 J for all layers. The tails of distributions look 
more evident due the log-scale of the vertical axis. 



The comparison of turbulence intensities J for each of the six layers occasionally 
shows detectable differences. The distributions of the differences in SJ are given in 
Fig. [21 Numerical values for the distributions are listed in Table [T] Relative fractions 
of 5 J fallen into left and right wings are presented in the two last rows. 



Table 1 Layer-by-layer comparison of the profile restoration using P and Nl techniques. 
Medians are given in units of 10~^'*m^'^ 



Layer altitude, 


km 




0.5 


1 


2 


4 


8 


16 


Median for Nl 






1.64 


1.33 


0.00 


2.43 


2.37 


2.36 


Median for P 






2.52 


0.52 


0.20 


2.91 


1.89 


2.42 


5J < -2.0- 10- 


-14 


% 


13.7 


4.2 


7.8 


10.9 


0.8 





5J > 4-2.0- 10- 


-14 


% 


8.2 


19.5 


6.0 


4.0 


6.6 






For low layers (0.5, 1 and 2 km) the base of their distribution is clearly asymmet- 
rical. This leads to the redistribution of the turbulence power between the adjacent 
layers and changes the median of layer intensity. In the 4 km and 8 km layers the effect 
is fainter and insignificant for the 16 km layer. It is apparent that the redistribution 
of the turbulence energy is related to the termination of the minimization process into 
certain local (not global) minimum, likely into the nearest starting point. 

This assumption is confirmed by Fig. [S] where the residuals of all the processed 
data are shown. Since the systems of equations are exactly the same, the cases where 
Rp > Rni (the points below the diagonal line) are occurrences of the algorithm P 




Fig. 3 Comparison of the solution residuals R^ for both methods (ai 50 000 points) 



finishing in a local minimum. For instance, in 50% of cases Rp is greater than 1.06-7?^]^ . 
In a few cases, the direct minimization doesn't converge in a finite number of iterations, 
and the profiles of such points are therefore absent. The number of such cases is less 
than 0.2%. 



5 Analysis of residuals R 



Analysis of residuals usually helps to determine the agreement between real data and 
the model. Since the R is the weighed sum of the squares of residuals, the direct 
comparison with the % distribution is possible. The distribution of values R is given 
on the left-hand plot of Fig. 3] The theoretical x distribution with 4 degrees of freedom 
(10 equations with 6 variables) is also plotted. One can see that the distribution in the 
case of Nl is closer to the theoretical curve but a significant part of solutions do not 
correspond to x criteria and should accordingly be dropped. In theory, x > 13.3 
must occur only in 1% of events, but really it is in 10% for Nl and 13% for P. That 
is why a large threshold of order 100 was used in the output filtering of the results by 
values of R^. The filter R^ > 100 rejects 2.2% for Nl and 3.8% in the case of P. 

Additionally, to show the adequacy of the model using R , the following conditions 
must be fulfilled: 1) the errors of the input data (vector b) must be estimated correctly, 
2) these errors are distributed normally and almost independent, 3) the number of 
degrees of freedom is determined correctly. 

First, the degrees of freedom are not known in advance for the linear problem with 
restrictions. In a sense, the restrictions are the equivalent of adding new equations to 
the system. The number of degrees of freedom increases when the restrictions come 
into effect (the solution is on the boundary of allowed area). The matrix rank on the 
area border is less than the rank of the source matrix. The NNLS algorithm excludes 
the matrix column which corresponds to a layers with zero intensity ("empty" layer). 




Fig. 4 Left: normalized distribution of R^ for the Nl (black) and P algorithms (grey line). 
The thin line is the x^ distribution with 4 degrees of freedom . Right: Relationship between 
median R^ and the total turbulence Jtotal in case of algorithm P (grey points), Nl (black 
points) and N2 (open circles), thin line denotes the mean value of degrees of freedom m 



The dependence on the total turbulence intensity Jtotal of the R medians and the 
mean of degrees of freedom m (number of empty layers+4) are plotted to emphasize 
its importance. Medians and means were calculated using bins of 1001 points of the 
Jtotal array. The median value of Jtotal is 1.85 • 10~ m ' , and 90 percentile is 5.41 • 
10-l3ml/3. 

These dependencies are shown on the right-hand plot of Fig.[4l In the case of strong 
OT (when in each layer the intensity is much more that error of solution) m should be 
4. In reality, m is about 5. The number m increases to 7 when the turbulence is weak, 
due to the appearance of "empty" layers. If R were exclusively dependent on m, the 
R distribution would match the m distribution shifted down by « 2/3. 

There must therefore be additional factors that decrease R with strong turbulence 
and increase it in case of weak one. The first such factor is the adequacy of the esti- 
mation of errors of input data and consequently the correctness of weighting the linear 
system. Obviously, that overshooting of errors leads to a decrease of R and vice versa. 

The systematic errors of the scintillation indices connected with the wrong MASS 
instrumental parameters or incorrect sky background subtraction or fight scattering 
[6] can be present in the input data. The influence of such errors as well as errors of 
weight functions Q{h) (inaccuracy of matrix A) is estimated in a different way. 



6 Properties of random errors of the input data 



The input data for the vertical profile restoration are 10 average scintillation indices 
(s ) calculated over an integration time Ta from individual measurements s . The error 
estimates of {s ) used for calculation of diagonal weight matrix W are computed in 
the ordinary way. The application of such weighting is valid if the input data errors 
are uncorrelated. 

This assumption is quite valid if the input errors are originated from the photon 
statistics and the scintillation process itself is homogeneous on the time scale of Ta- To 
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Table 2 Median values of variances of errors 


of mean scintillation indices 






Indices 


4 


2 


„2 „2 


2 
'^AB 


„2 


^AD 


<i2 


«2 


2 


<72 X 106 


5.26 


2.42 


0.86 0.17 


0.50 


1.16 


3.47 


0.19 


1.25 


0.20 



check these predicates, variances a^2 and correlation coefficients p^ were calculated 
over sets of indices s over Ta interval. Over a small sample size (n — Ta/Tf, = 60), 
the relative accuracy of the estimates is not better than 20%. The correlation and 
non-normality of the scintillation indices even more degrade the accuracy. 

The median errors o-,2\ of mean indices defined as cr ^ijn are given in Table [J] One 
can see that the maximal weights will be attributed to equations for indices D, BC and 
CD, and minimal weights for A and AD ones. 

Examination of the correlation coefficients pij shows that 11 of 45 median coef- 
ficients are greater than 0.8 and 21 are greater than 0.6. The correlations AB;C and 
AB:D are near zero which is evidence of the domination of random noise in scintilla- 
tion index s^^ and of weak connection via turbulent layers. On the another hand, the 
correlation between s^ and Sq is 0.88 which shows the small contribution of random 
noise in the indices. 

It is clear that those median values correspond to situations when Jtotal is close to 
its median value 1.85T0~ rn ' . Both the variances and correlation coefficients depend 
on the turbulence intensity. As expected, the values Py decrease if the turbulence calms 
down, and hence the fraction of random fiuctuation in indices is increased. 

Hence, it is the turbulence non-stationarity during the time Ta and not the photon 
noise introduces the valuable part in input data errors. It is easy to explain because 
during 1 minute the atmosphere drifts by 1 km with the wind velocity is about 20 m/s 
what is much more than the turbulence outer scale. 

It is known that, in the case of correlated input data errors, the best unbiased 
estimator of solution of linear system Q is the solution of equations weighted in the 
following manner: 

B~^Ax = B"^b, (5) 

where B is the triangle matrix obtained from the decomposition of the symmetric, 
non-negatively defined covariance matrix E of errors of the input data b. Such de- 
composition E = BB can always be performed, for instance by the Cholesky method 

m- 

The estimator of covariance matrix S is calculated from the same sample of indices 
that is used for determination of means (s ). As before, it leads to large errors in the 
weight matrix B~ . The uncertainty of the weights affects the solution weakly, but 
it does give a valuable estimation of the solution errors. The Cholesky decomposition 
cannot be performed in some rare cases. A damping coefficient of 0.8 is used for non- 
diagonal matrix E elements in order to prevent such cases. 

Having been weighted according ((5| the right hand side of the system holds uncorre- 
lated errors. This kind of modification (decorrelation of input errors) was implemented 
in the version N2 of algorithm. 



7 Estimation of random errors of the solution 

An estimation of the random errors of the solution can be obtained with a general 
method based on Gauss-Markov theorem [10) by calculating the unsealed covariance 
matrix: 

C = (ArA)-^ (6) 

The variance of k-ih layer error is defined as follow 

where R is the residual of the solution of correctly weighted set of equations, k is 
a number of equations (10 independent indices), r is the solution rank (a number of 
nonzero layers). The rows and columns of the matrix C corresponding to "empty" 
layers are identical to zero. 

An estimation of correlation coefficients pij are obtained directly from the covari- 
ance matrix 

Pij ^ij / [^ii^jj ) • \^} 

The above relationships assume the correct covariance matrix S, that is satisfied 
only in algorithm N2. Nevertheless, the calculation of solution errors was implemented 
in algorithm Nl in order to estimate an effect of applying formula ([5]). 



8 Solutions and their errors for algorithms Nl and N2 

The cumulative distributions of layer by layer OT obtained with algorithms Nl and 
N2 are given on the left-hand plot of Fig. [5] One can see that the curves are very close. 
The maximal difference is observed in 1 km layer but it doesn't exceed 3 • 10~ rn ' . 
Both curves are similar to the cumulative distribution in paper [7] which was obtained 
with algorithm P. 

The behavior of residuals R obtained for the correctly weighted system ((5| is given 
on the right-hand plot of Fig. [l] The curve is significantly higher and has a smaller 
inclination than one would arise from theory. The considerable increase of the residuals 
with decreasing Jfotal could mean that the systematic discrepancy (of additive type) 
between the input data and the model begins to dominate. However, as it will be 
shown, the systematic effects are significantly lower. 

For each layer, the dependence of median of the standard deviation ctj on Jfotal were 
studied in order to determine the typical errors of the restored profile. For algorithm 
Nl this dependence is shown on the right-hand plot of Fig. (5] As before, the median 
was calculated by the sample of 1001 points. The data related to empty layers were 
ignored. This method gives an estimate of the upper limit of the error. 

Notice that the layers are divided into three groups: 1) the 16 km layer is determined 
with precision better than 5% even in case of weak turbulence 2) layers 8 and 4 km are 
determined with precision better than 10% when Pfree > 0.25" and close to 10% in 
case of weaker turbulence 3) the layers 2, 1 and 0.5 km are determined with precision 
worse than 10% when /Jf^ee < 0.7" and worse than 20% when Pfj-ee < 0.3". 

The estimations of Oi axe slightly greater for algorithm N2 than for Nl especially 
in case of strong turbulence. However, they are close in case of weak turbulence. Ex- 
trapolated to zero turbulence (Jj are respectively 1.3, 0.8, 0.7, 0.3, 0.2, 0.1 in units of 
10-1^1/3. 
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Fig. 5 Left: cumulative distribution of the OT intensity in 6 layers for algorithm N2 (lines) 
and Nl (dashed lines). Right: dependence of medians of the solution errors on Jtotal for the 
Nl algorithm. The thin lines are the percentages of Jtotal 



For the algorithm Nl, the correlations of errors pij in the nonempty, nearest- 
neighbor, layers are approximately —0.9 to —0.95. For the algorithm N2, the values of 
p were expected to decrease significantly but, instead, only did so by a small amount. 
This means that the structure of the equation set and not the characteristics of the 
input data is the dominant source of error correlation. The correlation of errors de- 
creases with increasing distance between the layers. For instance, for the 2 km and 
16 km layers, p is between -1-0.5 and —0.5 depending on the number of empty layers 
between them. 

The sum of the intensities of the layers is very stable thanks to the error correlation. 
The effect of "dragging" of the turbulence to the nearest-neighbor layers was marked 
earlier [18]. Note that such a strong correlation is the consequence of solving, and not 
of the physics of the phenomenon. 

The question about the precision of the OT measurements in each layer was raised 
in [5] , in the form of estimation of the noise of the restoration process. Generated by the 
turbulence motion over intervals less than Ta errors are included in this noise and in 
agreement with the ones given here by order of magnitude and behavior. Our estimates 
coincide with the estimates of errors in [IB] as well. 

The problem with making valid estimates of the error in each layer is of great im- 
portance because the value of the error can be comparable with the intensity of the 
OT in the layer, especially for lower layers (see Fig. [5]). It leads to significant widen- 
ing of Ji distribution and the underestimation of both medians and minor percentiles. 
The requirement of non-negativeness leads to an increase of zeros in the cumulative 
distribution. In some cases the median turbulence in a layer becomes zero. The un- 
derestimation of the intensity of low-altitude turbulence can significantly change the 
results of AO simulation. Fig. [5] shows that only the 16 km layer has a clearly defined 
distribution. 

In principle, knowing the estimate for ai, one can try to reestablish the actual 
distribution of turbulence using deconvolution, but the problem is not trivial because 
the distribution of errors depends on the turbulence in other layers. 
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Fig. 6 Left: cumulative distributions of the turbulence intensity Ji in 6 layers after restoration 
witli 1 s intervals foUowred by a 1 min averaging. The distributions with a 2 min averaging are 
shown by the dashed lines. Right: dependencies of the median errors (Ji on the total turbulence 
Jtotal for N3 algorithm 



9 The final version: algorithm N3 



The significant influence of errors on the shape of the distribution rises the question: 
how does the distribution depend on the integration time Ta and what happens during 
the averaging of near-neighbor data? It is clear that, in the linear problem without 
restrictions, the averaging of the input data over Ta is equivalent to the averaging of 
output data over the same interval. 

Detailed research of temporal behavior of the atmospheric turbulence intensity 
is beyond the scope of this paper. Rough estimates show that the auto-correlation 
function of the scintillation process can be fitted with ~ exp (— i/io)i where to is about 
tens of minutes that is much more than the integration time. This way, measurements 
of turbulence that are close in time correlate strongly (for details, see the technique 
of noise estimation in [5]). This leads to the fact that changing the accumulation 
time during the profile restoration or the following averaging do not alter the results 
statistics, when the noise is negligible. The well-determined Jig can be given as an 
example. Its distribution practically doesn't depend on Ta (the difference is less than 
1%) between 5 and 240 s or by averaging 2, 4, 8 and 16 consecutive points. 

For the other layers where the noise is significant, the cumulative distributions 
change according to integration time and averaging, owing to the non-linearity of the 
problem under discussion. This was the reason to change the averaging sequence in 
algorithm N3. 

In the final version N3, the covariance matrix E is calculated over all the data 
within Ta but the OT profiles are restored using individual indices s (integration time 
Tfy). The obtained profiles are averaged over an accumulation interval Ta after that. 
Such algorithmic approach only becomes conceivable after increasing the speed of the 
computational processing. The cumulative distributions obtained with algorithm N3 
are given on the left-hand plot of Fig. [6l 

It can be seen that the distributions of all of layers (with the exception of the 16 km 
layer) have undergone significant changes relative to the curves in Fig. [5] A jump in 
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zero for 0.5 km layers disappeared almost entirely. The medians increase for all of the 
layers, nonzero 25% percentile appears for other layers. 

The median of the turbulence in layer 0.5 km is 5.95 ■ 10~ m ' , that is more 
that 3 times the former value. The overall turbulence Jtotal has increased slightly and 
its median value becomes 2.12 • 10~ rn ' , in other words it is 13% greater than the 
algorithm Nl gave, the seeing Pfree grew from 0.46" to 0.50". 

The question about the statistical adequacy of the output set is of great importance 
because the distribution of the individual solutions is extremely asymmetrical in case 
of weak turbulence due to the requirement of non-negativeness. Usage the median as 
such a characteristic leads to a turbulence distribution similar to the ones shown in 
Fig. O Moreover, in this case the sum of the layer intensities can differ greatly from 
the turbulence integrated over the whole free atmosphere. 

Clearly, the sum of medians may differ from median of all the turbulence. But if 
this difference is large, the appropriateness of such an estimator is in question. In the 
case of the algorithm Nl the sum of medians is only 55% of the total turbulence, while 
the algorithm N3 gives 75%. 

The usage of a mean as a central estimator leads to a physically more understand- 
able result: the probability of zero turbulence must be zero. That estimator doesn't have 
to be unbiased, because the noise leads to the appearance of non-zero mean value even 
in the layer with zero intensity. However, such kind of bias are intrinsic for any problem 
with restrictions. The mean is well matched to the median in case of well-determined 
turbulence and both of these estimators are unbiased in practice, for instance for high 
layers. 

Estimates of the error of the mean over an accumulation time Ta are calculated 
directly as root-mean-square. The behavior of the errors calculated in such a way is 
shown on the right-hand plot of Fig. |6l It can be seen that the errors decrease by a 
factor of 3 to 5 compared to the N2 algorithm. 

The accuracy of the sample variance is not high because of the amount of data in 
the 1 min integration, and the excess coefficient 72 for intensity distribution varies from 
for the layers with fully sustained turbulence to 10 in the case where they are almost 
"empty". The relationship for the relative error of the variance ££,.2 = v(72 + 2)/]V 
shows that the accuracy of the estimator varies between 15 and 40 %. While calculating 
the sample variance, the intensities are also assumed to be uncorrelated on account of 
their noticeable noise. 

Note that in the case of N3, the criterion of the restoration quality is the average 
value of R of individual solutions, which is distributed almost normally. Hence, the 
traditional way of detecting non- valid results is possible. 



10 Possible methods for turbulence profile discretization 

Let hi be the altitude grid and J^ the turbulence intensity collected @ between layers 
hi and /i^+i. Strictly speaking, we obtain the following equation from the mean- value 
theorem: 

As^ = / '^' Cl(h)Q{h)dh = J.,Q{hi), (9) 

Jhi 

where hi is some altitude between hi and /i^+i. Moreover, it will differ for the different 
indices Sj. The equation set ([3]) assumes hi j are the same for all indices and matches 
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Fig. 7 The restoration process; contribution of tiie tiiin turbulent layers of unitary intensity at 
altitude h. The thick dashed curves are the solution residuals VR?, the solid line is a sum over 
all layers. The thin dashed curves are intensities in the corresponding layer. Left: point-like 
representation. Right; with constant C^ in the layer 



some effective altitude h* to which Ji is bound; for instance, a mean or geometrical 
mean. For this model the real layer thickness is not known exactly which doesn't lead 
to a problem when operating in Jj space, but estimation of the corresponding C„(ft) 
becomes inaccurate. 

Little error is introduced if Cn{h) and Q{h) do not change rapidly inside the layer, 
i.e. the layer is thin: Ah/h <C 1. The correspondence hi = h* is also valid in the 
case where all the turbulence is located at the effective altitude. In reality, the layer 
thickness is comparable with its altitude, because hi ~ 2/ij_i is necessary to fill all the 
altitude range with 6 layers. Nevertheless both Q{h) and Cn{h) can change a number 
of times within the layer. 

The second trivial representation is based on the assumption that the turbulence 



Cn(h) is constant within the layer and Ji — Cn{hi){h. 



i+l 



hi). Then 



As^ = 



hi 



Cn{h)Q{h)dh ^ 



■h 



hi 



hi 



hi 



Q{h)dh = J,{Q{h)), 



(10) 



and the matrix of the equation set has the mean values of the weighting function in 
layer {Q{h)) instead of the values Q{h*). 

The third possible form is piecewise linear continuous function for the structural 
coefficient: Cn{h){hi^i — hi) = Cn{hi){hi^i — h) + Cn{hi^i){h — hi). The turbulence 
Cn{hpf) = at the top border of atmosphere (~ 25 -^ 30 km). The coefficients of 
equation set are also calculated from the integrals of the functions Q{h) and hQ{h) 
over the layer. Unlike previous representations, a composing the system for Cn{hi) in 
grid nodes is preferred. 

The contribution of each turbulence layer with unitary intensity located on altitude 
h was computed for all of the described representations on the grid with nodes at 0.5, 
1, 2, 4, 8 and 16 km. Similar contribution was considered in [5], but the goal of that 
research was to show the qualitative behavior and the possibility of the solution but 
not the examination from a mathematical viewpoint. The results of modelirrg of the 
first two paths are given on Fig. [T] 

The right-hand part of the equation set corresponds to the values of the weighting 
functions Q{h) at altitude h without noise. In this case the problem is kept linear and 
can be solved without restriction, nevertheless, the algorithm Nl was applied (N2 and 
N3 were developed at a later stage). 
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Fig. 8 Top: Reaction of the N3 algorithm to a thin turbulent layer of intensity J = 5 • 
10~^^ m^' ^ at altitude h and real errors of the scintillation indices from night of September 
17, 2005. The dashed line is the expected total intensity, the thin black is the measured total 
intensity. Bottom: behavior of the residual R^. Left: 6-layer restoration model. Right: 12-layer 
restoration model 



A characteristic of the solution's residuals is evident: in the first version, the residual 
becomes zero when the layer altitude coincides with a node altitude. For the second 
and third representations, the residual is never and the altitudes of minimums do not 
match the nodes exactly. 

The behavior of the residual changes mainly when the test layer is higher than the 
upper node (16 km in our case). It increases significantly and reaches 2.5 at 20 km 
(becomes more than 6 in sense of _R ) . The main conclusion from this fact is that the 
model isn't adequate in case of existence of noticeable turbulence at altitudes above 
16 km. 

The restored total turbulence intensity J exceeds unity by 0.05 hence the algorithm 
overestimate the total turbulence. All the turbulence above the upper node is assigned 
to it in the intensified form. For instance the turbulence at 25 km will be appended to 
16 km node with a factor of 1.5. Significant turbulence at up to 25 km was noted in 
a number of researches. For instance, the SCIDAR measurements [19] show the mean 
summer intensity above 16 km to exceed 1.5 ■ 10^ m ' . 

It is possible to combine the requirements of increasing the upper boundary and 
decreasing the layer thickness by transitioning from a 6-layer model to models with a 
larger number of layers. For example, by adding a lower node at an altitude close to 
0.35 km and an upper node at an altitude of about 25 km in order to guarantee the 
absence of turbulence above the last node. 

The differences between the three representations of the problem are not important 
and there is no need to reject the linear system of the form ^. 

In the case where noise is present, the reaction of the model to thin turbulent 
layers is different. In order to simulate a real situation, the indices s were normaUzed 
in such a way that their average over Ta corresponds to a turbulence layer with fixed 
intensity at a needed altitude. The relative indices fluctuations are kept in this case. 
If the photon noise does not dominate then the covariance matrix and the weighting 
matrix are unchanged. 
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Fig. 9 Processing results of the measurements made at Maidanak on September 17 2005. 
Left: 6-layer restoration model. Right: 12-layer model. Bottom panels: free seeing fifree 



Simulation results are given on the left-hand plot of Fig. [S] The measurements of 
September 17, 2005 were used as noise template. One can see two significant differences 
from the noiseless model: 1) The residual is much closer to its average value if the layer 
is located below the 0.5 km node (the first node of the altitude grid) 2) The residual 
R increases rapidly with increasing altitude of the layer above the upper node, but 
the solution continues to fall. Similar behaviors were observed for the entire night of 
data. 

One can see that the mean value of Jtotal is about 10% over the input intensity 
partially due to the noise bias. Also, the intensity in the first and the last layers is 
slightly overestimated. The "tail" of the 0.5 km layer originated by upper layers noise 
is a characteristic detail. A dim star (~ 5.5 pulse/ms) was acquired in the second half 
of the night and the intensity of the noise "tail" increased to about 30% of Jtotal- 
However, if a layer is added at 0.35 km , all the noise migrates to it and the 0.5 km 
curve becomes normal. 



11 Study of the profile restoration Avith larger number of layers 



Efforts to obtain OT profiles with more than 6 layers were previously attempted with 
the algorithm P, but even the addition of two layers resulted in a significant growth 
of the solution's instability. Clearly, in the case of an unrestricted linear problem, 
the maximal size of the solution cannot be larger than the number of input data (10 
scintillation indices in our case). The restrictions give the possibility to increase the 
number of grid nodes above 10. The new algorithm N3 which uses the great regularizing 
properties of NNLS makes it possible to use twice the number of layers without loss of 
solution accuracy. 

Note that the denser grid does not necessarily increase of resolution. For a real 
resolution increase, a low level of input "atmospheric" noise (stationarity of OT) is 
needed too. In any case, additional nodes let to localize features of the turbulence in 
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Fig. 10 Left: Median turbulence intensity profiles at Maidanak for both 6-layer (open circle) 
and 12-layer (black circles) models. Asterisk — 12-layer data recalculated on 6 layers for 
comparison, grey crosses — data from 7 . Right: Median turbulence intensity profiles at Tolar 
in January 2004 (black points) and in October 2005 (open circles) obtained with 12-layers 
model. Dashed line — the mean profile 



altitude more precisely. It is very important for understanding a turbulence generation 
process above a site. The study of a 12-layer restoration is given below although the 
experiment shows that solving with 14-15 layers is possible too. 

For simplicity the 12-layer grid was obtained from the 6-layer grid by adding in- 
termediate layers with 1.5 times the altitude. In Fig. [^ the results of 6 and 12-layers 
OT restoration are shown. The night of September 17 2005 is calm enough, although 
a number of bursts are registered during the night practically in each layer. The com- 
parison shows that high turbulence is located not at 16 km but between 8 and 12 km 
where the tropopause must be. 

The total integrated turbulence is the same as one can see from comparison of see- 
ings Pfree Computed by profiles. In the 22 km layer the weak turbulence is collected, 
its reality may be questioned but this layer is needed to accumulate the highest tur- 
bulence. Low altitude turbulence is attributed to the 0.5 km layer. When 13-th layer 
at altitude 0.35 km was added, majority of 0.5 km turbulence migrates to this layer. 
Most likely, the turbulence observed at 0.5 km comes mainly from lower altitudes. 

Turbulence intensity profiles for 6 and 12-layers models are shown in Fig. 1101 It 
can be seen that they have a physically correct behavior: strong near the ground and 
boundary turbulence become apparent in 0.5 km layer, next increase at 6 - 10 km in 
tropopause zone and decrease at highest altitude after. 

To estimate the real altitudinal resolution, the thin layer generation described above 
was used. As for the 6-layers model the night September 17 2005 was selected as a 
noise template. The simulation results are presented on the right-hand plot of Fig. [H] 
In general, the response looks similarly to the 6-layer model, however the lowest and 
the highest layers differ more from the others. The maximal intensity in intermediate 
layers is « 75% of the theoretical value. This indicates that the 12-layer grid is denser 
than the resolution for the current "atmospheric" noise produced by non-stationarity. 

The picture is more complex for real turbulence distribution where a strong layer 
may be near weak one. Depending on their disposition, the weak layer turbulence flows 
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Fig. 11 Results of the OT profiles restoration for data of measurements on Cerro Tolar on 
October 1 (left) and October 19 (right) 2005 



partially either to or from the strong layer. Of course, the 12-layer model is more 
sensitive to this effect and to the noise than the 6-layer model. 

The problem of the top layer is solved by adding an extra node at an altitude of 23 
- 25 km, where the OT is very weak. Dealing with the bottom layer is a more complex 
problem. In this layer the turbulence is very strong and it therefore accumulates a lot 
of noise. An additional node makes it possible to specify the lower boundary of the 
OT. For example, in Fig. [TT] the layer at 0.35 km drives the turbulence higher than 
0.4 km to the 0.5 km layer which is free from most fraction of the noise power. 



12 Validation of the profile restoration for MASS/DIMM device 
measurements 



Previous studies were performed using the original MASS instrument. However, most 
of the data are delivered nowadays by MASS/DIMM devices [2j, which have differ- 
ent aperture set and spectral response — i. e. different set of atmospheric weighting 
functions Qj{h). 

Although these differences are not fundamental, it was necessary to verify the 
final algorithm N3 for the MASS/DIMM data. The measurements at mount Cerro 
Tolar (Chile) carried out by the TMT group ^KWIT\ in January 2004 (11 nights) and 
October 2005 (17 nights) were used as input data. In January, the MASS/DIMM 
segmentator mask produced a ghost image in the B channel direction, which resulted 
in large residuals in the restoration. The problem was fully fixed by accounting for 
scattered light as large as 0.12. 
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Profiles were restored using 12-layers grid. Tiie character of tlie OT in tfiese two 
seasons was found essentially different (see right-hand plot on Fig. I10|) . In Fig. Illl the 
restored OT in October 1 and October 19 nights are shown. 

To detect the accumulation of "undue" turbulence in the lowest and highest layers, 
two very low layers (0.25 and 0.35 km) were added to the 12-layers (0.5, 0.7, 1 ... 16, 22 km) 
model. Additionally, in the processing of October 19 a very high layer 32 km was added. 

One can see that the N3 restoration algorithm works perfectly for both the 14 and 
15 layer grid. On October 1, residual R is about 1 before 4 UT and « 2 after. The 
reason is that a bright star (750 pulses/ms in D aperture) was changed by a fainter 
one (« 200 pulses/ms) in that moment. This event clearly affects the 0.25 km layer 
intensities. On October 19 a bright star was used until 2 45™ UT, therefore the strong 
0.25 km turbulence before 3 15™ is certainly real. 

Measurements made with the SLODAR 2_2 show that stratification is often ob- 
served in the boundary layer. So, weaker turbulence in the 0.35 km than in its neigh- 
boring layers may be real. Also, on can see that uppermost OT is not located at 16 km 
but in the 8-12 km altitude range. Such correction increases the isoplanatic angle 
when calculated from the profiles. 



13 Conclusion 

The described modifications of the OT profiles restoration process do not, in principle, 
change the physical results obtained with MASS instrument. However, the revised 
restoration algorithm makes the results more reliable and clearer for interpretation. 
Studies of the profile restoration were carried out before in connection with MASS 
instrument verification [5ll6lll8) . Many discovered effects stimulated the work presented 
in this paper. 

We formulate a summary of the modifications made to the algorithm and their 
consequences on the results obtained with the MASS technique. 

1. Algorithm based on direct minimization is replaced by a more mathematically 
correct method of Non-Negative Least Squares (NNLS). 

2. Comparative analysis of the restoration of altitude profiles of the optical turbu- 
lence on measurements at mount Maidanak in 2005 - 2007 shows that the NNLS 
algorithm works always better. 

3. Calculation of solution errors has shown that in the case of weak turbulence errors 
are comparable with the turbulence intensities. 

4. The consequence of that is the distortion of the statistical distribution of layer 
intensities, leading to a significant underestimation of the median and quartile 
values. 

5. The main source of uncertainty is the non-stationarity of OT at the scale of the 
integration time (usually 60 s), noise of the measurements is almost always much 
lower. 

6. De-correlation of the input data errors slightly changes the estimates of output 
errors and does not change the solutions themselves. 

7. The change of the strategy of the restoration algorithm from "averaging indices - 
restoration" to "restoration - averaging profiles" has led to a significant decrease 
of the errors of the OT profiles. 
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8. As a result, vertical profiles have more statistical significance and physical clarity. 
The integrated turbulence intensity in free atmosphere has increased by no more 
than 15 %. 

9. In the case of turbulence, corresponding to good seeing, it is possible to increase 
the vertical resolution and accuracy of turbulent layer localization (from h/2 to 
h/4) using the 12 - 14 layers model. 

10. The new algorithm is also effective for data obtained with the MASS/DIMM in- 
strument, which has different set of the entrance apertures. 

An additional consequence from the revision of restoration process is the improve- 
ment of the C-|— I- code that is clearer and more suitable for further modifications and 
a substantial (order of magnitude) increase in processing speed. 
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